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Abstract. This paper studies mixed finite element approximations of the viscosity solution to 
the Dirichlet problem for the fully nonlinear Monge-Ampere equation det(D 2 M°) = / (> 0) based 
on the vanishing moment method which was proposed recently by the authors in [19| . In this 
approach, the second order fully nonlinear Monge-Ampere equation is approximated by the fourth 
order quasilinear equation — eA 2 « £ + det D 2 u e = f. It was proved in [17[ that the solution u e 
converges to the unique convex viscosity solution u° of the Dirichlet problem for the Monge-Ampere 
equation. This result then opens a door for constructing convergent finite clement methods for the 
fully nonlinear second order equations, a task which has been impracticable before. The goal of 
this paper is threefold. First, we develop a family of Hermann-Miyoshi type mixed finite element 
methods for approximating the solution u E of the regularized fourth order problem, which computes 
simultaneously u e and the moment tensor a e := D 2 u e . Second, we derive error estimates, which 
track explicitly the dependence of the error constants on the parameter e, for the errors u e — ui 
and <r e — cr|. Finally, we present a detailed numerical study on the rates of convergence in terms of 
powers of e for the error u° — u E h and o e — o e h , and numerically examine what is the "best" mesh size 
h in relation to e in order to achieve these rates. Due to the strong nonlinearity of the underlying 
equation, the standard perturbation argument for error analysis of finite element approximations of 
nonlinear problems does not work for the problem. To overcome the difficulty, we employ a fixed 
point technique which strongly relies on the stability of the linearized problem and its mixed finite 
clement approximations. 
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1. Introduction. This paper is the second in a sequence (cf. [20]) which con- 
cerns with finite element approximations of viscosity solutions of the following Dirich- 
let problem for the fully nonlinear Monge-Ampere equation (cf. [2"3"]1: 



where is a convex domain with smooth boundary SO. D 2 u°(x) and det(D 2 u°(x)) 
denote the Hessian of u° at x £ and the determinant of D 2 vP(x). 

The Monge-Ampere equation is a prototype of fully nonlinear second order PDEs 
which have a general form 



with F(D 2 u° , Du° , u° , x) — det(D 2 u°) — f. The Monge-Ampere equation arises nat- 
urally from differential geometry and from applications such as mass transportation, 
meteorology, and geostrophic fluid dynamics [U [8]. It is well-known that for non- 
strictly convex domain f2 the above problem does not have classical solutions in gen- 
eral even /, g and d£l are smooth (see 22J). Classical result of A. D. Aleksandrov 
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(1.1) 

(1.2) 



det(L> 2 u°) = / 



in il C R™, 
on d£l, 




(1.3) 



F(D 2 u Q ,Du a ,u°,x) = 
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states that the Dirichlct problem with / > has a unique generalized solution in the 



class of convex functions (cf. [SIS])- Major progress on analysis of problem (1.1 1 



(1.2) has been made later after the introduction and establishment of the viscosity 
solution theory (cf. US])- We recall that the notion of viscosity solutions was 
first introduced by Crandall and Lions |llj in 1983 for the first order fully nonlinear 
Hamilton-Jacobi equations. It was quickly extended to second order fully nonlinear 
PDEs, with dramatic consequences in the wake of a breakthrough of Jensen's max- 
imum principle [25] and the Ishii's discovery [53] that the classical Perron's method 
could be used to infer existence of viscosity solutions. To continue our discussion, 
we need to recall the definition of viscosity solutions for the Dirichlct Monge- Ampere 



problem (|l.l|-(|1.2| (cf. [23]). 

Definition 1.1. a convex function u° £ C°(£!) satisfying u° = g on dfl is called 
a viscosity subsolution (resp. viscosity supersolution) of ( |1.1[ ) if for any ip £ C 2 there 
holds det(D 2 (p(xo)) < f(xo) (resp. det(_D 2 ip(xo)) > f(xo)) provided that u° — ip has 
a local maximum (resp. a local minimum) at xq G f2. u° E C°(fl) is called a viscosity 
solution if it is both a viscosity subsolution and a viscosity supersolution. 

It is clear that the notion of viscosity solutions is not variational. It is based 
on a 11 differentiation by parts 11 approach, instead of the more familiar integration 
by parts approach. As a result, it is not possible to directly approximate viscosity 
solutions using Galerkin type numerical methods such as finite element, spectral and 
discontinuous Galerkin methods, which all are based on variational formulations of 
PDEs. The situation also presents a big challenge and paradox for the numerical PDE 
community, since, on one hand, the 11 differentiation by parts" approach has worked 
remarkably well for establishing the viscosity solution theory for fully nonlinear second 
order PDEs in the past two decades; on the other hand, it is extremely difficult (if all 
possible) to mimic this approach at the discrete level. It should be noted that unlike 
in the case of fully nonlinear first order PDEs, the terminology "viscosity solution" 
loses its original meaning in the case of fully nonlinear second order PDEs. 

Motivated by this difficulty and by the goal of developing convergent Galerkin type 
numerical methods for fully nonlinear second order PDEs, very recently we proposed 
in |17] a new notion of weak solutions, called moment solutions, which is defined 
using a constructive method, called the vanishing moment method. The main idea 
of the vanishing moment method is to approximate a fully nonlinear second order 
PDE by a quasilinear higher order PDE. The notion of moment solutions and the 
vanishing moment method are natural generalizations of the original definition of 
viscosity solutions and the vanishing viscosity method introduced for the Hamilton- 
Jacobi equations in [11]. We now briefly recall the definitions of moment solutions 
and the vanishing moment method, and refer the reader to |17[ 119] for a detailed 
exposition. 

The first step of the vanishing moment method is to approximate the fully non- 



linear equation (1.3) by the following quasilinear fourth order PDE: 

(1.4) -sA 2 u £ + F(D 2 u e ,Du £ ,u £ ,x) = (e > 0), 

which holds in domain fl. Suppose the Dirichlct boundary condition m° = g is pre- 
scribed on the boundary d£l, then it is natural to impose the same boundary condition 
on u e , that is, 

(1.5) u £ — g on 9f2. 



However, boundary condition (1.5 1 alone is not sufficient to ensure uniqueness for 



fourth order PDEs. An additional boundary condition must be imposed. In |17] the 
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authors proposed to use one of the following (extra) boundary conditions: 
(1.6) Au £ = e, or D 2 u £ v -v = e on dtt, 

where v stands for the unit outward normal to dfl. Although both boundary condi- 
tions work well numerically the first boundary condition Au £ = e is more convenient 
for standard finite element methods, spectral and discontinuous Galerkin methods (cf. 
[20]). while the second boundary condition D 2 u £ v ■ v = e fits better for mixed finite 
element methods, and hence, it will be used in this paper. 

In summary, the vanishing moment method involves approximating second or- 



der boundary value problem (1.2 1— (1.31 by fourth order boundary value problem 



(1.4 1— (1.5 1, (1.6 1. In the case of the Monge- Ampere equation, this means that we 
approximate boundary value problem ( |l.l[ )-( 1.2 1 by the following problem: 



(1.7) 
(1.8) 
(1.9) 



-eAV + det(D 2 u £ ) = f 
u £ = g 
D 2 u £ v ■ v = e 



in f2, 
on <9f2, 
on <9f2. 



It was proved in [T5] that if / > in f2 then problem ( |1.7[ )-( 1.9 ) has a unique solution 
u £ which is a strictly convex function over f2. Moreover, u £ uniformly converges 
as e — > to the unique viscosity solution of ( |l.l[ )-( [l~2| . As a result, this shows 



that (1.1 1 — ( 1.2 1 possesses a unique moment solution that coincides with the unique 



viscosity solution. Furthermore, it was proved that there hold the following a priori 
bounds which will be used frequently later in this paper: 



(1.10) 
(1.11) 



\ u Whj = 
\D 2 u £ \\ L2 



= 0(e 



\\U Hvt/2,00 = 
|| C Ol(Z>V)| 



0{e 



for j = 2,3. Where cof(D 2 u £ ) denotes the cofactor matrix of the Hessian, D 2 u £ . 

With the help of the vanishing moment methodology, the original difficult task of 
computing the unique convex viscosity solution of the fully nonlinear Monge-Ampere 



problem ( 1.1 )-( 1.2 ), which has multiple solutions (i.e. there arc non-convex solutions), 
is now reduced to a feasible task of computing the unique regular solution of the 



quasilinear fourth order problem ( 1.7 H 1-9 1. This then opens a door to let one use 



and/or adapt the wealthy amount of existing numerical methods, in particular, finite 
clement Galerkin methods to solve the original problem ( 1.1 1— ( 1.2 1 via the problem 

The goal of this paper is to construct and analyze a class of Hermann-Miyoshi 



type mixed finite element methods for approximating the solution of (1.7)— (1.9|. In 



particular, we are interested in deriving error bounds that exhibit explicit dependence 
on e. We note that finite element approximations of fourth order PDEs, in particular, 
the biharmonic equation, were carried out extensively in 1970's in the two-dimensional 
case (see [TU] and the references therein), and have attracted renewed interests lately 
for generalizing the well-know 2-D finite elements to the 3-D case (cf. [351 [3EJ I3~i] ) 
and for developing discontinuous Galerkin methods in all dimensions (cf . [TSl [57] ) . 
Clearly, all these methods can be readily adapted to discretize problem ( |1.7| -( [l~9| 
although their convergence analysis do not come easy due to the strong nonlinearity of 
the PDE (1.7 1. We refer the reader to [20, 28 for further discussions in this direction. 



A few attempts and results on numerical approximations of the Monge-Ampere 
as well as related equations have recently been reported in the literature. Oliker 
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and Prussner |30j constructed a finite difference scheme for computing Aleksandrov 
measure induced by D 2 u in 2-D and obtained the solution u of problem ( 1.7 1 — ( 1.9 1 



as a by-product. Baginski and Whitaker [5] proposed a finite difference scheme for 
Gauss curvature equation (cf. [19 and the references therein) in 2-D by mimicking 
the unique continuation method (used to prove existence of the PDE) at the discrete 
level. In a series of papers (cf. [T3| and the references therein) Dean and Glowinski 
proposed an augmented Lagrange multiplier method and a least squares method for 



problem (1.7 1— (1.91 and the Pucci's equation (cf. [Tj [55]) m 2-D by treating the 
Mongc- Ampere equation and Pucci's equation as a constraint and using a variational 
criterion to select a particular solution. Very recently, Obcrman 29J constructed some 
wide stencil finite difference scheme which fulfill the convergence criterion established 
by Barles and Souganidis in [3J for finite difference approximations of fully nonlinear 
second order PDEs. Consequently, the convergence of the proposed wide stencil finite 
difference scheme immediately follows from the general convergence framework of [3J. 
Numerical experiments results were reported in [SUJ [251 [13] , however, convergence 
analysis was not addressed except in |29j . 

The remainder of this paper is organized as follows. In Section [2) w e first de- 



rive the Hermann-Miyoshi mixed weak formulation for problem (1.7l-(1.9l and then 
present our mixed finite element methods based on this weak formulation. Section 
[3] is devoted to studying the linearization of problem ( 1.7 1- ( 1.9 1 and its mixed finite 



element approximations. The results of this section, which are of independent inter- 
ests in themselves, will play a crucial role in our error analysis for the mixed finite 
element introduced in Section [2] In Section [4] we establish error estimates in the 
energy norm for the proposed mixed finite clement methods. Our main ideas are to 
use a fixed point technique and to make strong use of the stability property of the 
linearized problem and its finite clement approximations, which all are established in 
Section |3] In addition, we derive the optimal order error estimate in the i/^-norm for 
u £ — u £ h using a duality argument. Finally, in Section [5] we first run some numeri- 
cal tests to validate our theoretical error estimate results, we then present a detailed 
computational study for determining the "best" choice of mesh size h in terms of s 
in order to achieve the optimal rates of convergence, and for estimating the rates of 
convergence for both u° — uf and u° — u e in terms of powers of e. 

We conclude this section by remarking that standard space notations are adopted 
in this paper, we refer to [SJ|22[TU] for their exact definitions. In addition, denotes 
a bounded domain in R™ for n = 2,3. (•, •) and (•, •) denote the L 2 -inner products on 
f2 and on <9f2, respectively. For a Banach space B, its dual space is denoted by B* . 
C is used to denote a generic e-independent positive constant. 



2. Formulation of mixed finite element methods. There are several popu- 
lar mixed formulations for fourth order problems (cf. [51 [TOJ US])- However, since the 
Hessian matrix, D 2 u e appears in ( 1.7 1 in a nonlinear fashion, we cannot use Au £ alone 
as our additional variables, but rather we are forced to use a e :— D 2 u e as a new vari- 
able. Because of this, we rule out the family of Ciarlct-Raviart mixed finite elements 
(cf. [10 ). On the other hand, this observation suggests to try Hermann-Miyoshi or 
Hermann- Johnson mixed elements (cf. jSJUS]), which both seek a e as an additional 
unknown. In this paper, we shall only focus on developing Hermann-Miyoshi type 
mixed methods. 
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We begin with a few more space notation: 

V := H 1 ^), W := {fi e [H 1 ^)}™"; fi l3 = 

V := H^Sl), V g := {v e V; w|an = 5}, 

W e := {(1 eW; fin- n\ dn = e}, W := € W; pin ■ n\ dn = 0} 



To define the Hermann-Miyoshi mixed formulation for problem (1.7 1- (1.9 1, we 
rewrite the PDE into the following system of second order equations: 



(2.1) 

(2.2) 



cr e - D 2 u £ = 0, 
-eAtr(a £ ) + det(<j £ ) = /. 



Testing (2.2 1 with v G V yields 

e I div(cr e ) • Dv dx + J det(a e )v dx — J fvdx. 
Jn Jn Jn 



(2.3) 



Multiplying (2.1) by fi G Wq and integrating over Q we get 



(2.4) 



f <7 £ :udx+ f Du £ ■ div(/i) dx — f fin-Tkjr^ds, 
Jn Jn fc=1 ./an c* 1 "*: 



where <7 £ : /1 denotes the matrix inner product and {r±(x), tq,(x) • • ■ , r n _i(a;)} denotes 
the standard basis for the tangent space to Oil at x. 

From ( 2.3 ) and ( 2.4 1 we define the variational formulation for ( 2. 1 )-( 2.2 ) as follows: 
Find (u £ , a e ) S V g x W E such that 



(2.5) 
(2.6) 

where 



(a £ , a) + (div( M ) ,Du E )=(g, fi) V/x G Wo, 
(div(a £ ),i;) + -(deta e , W ) = (/ e , V ) Vu G Vo, 



<5> M> = 



n—l 



fin ■ n 



and r = -f. 

e 



Remark 2.1. We note that det(a £ ) = ±§ £ D 2 u £ = 1 V™ . • < , for 7 = 
1,2, ...,7i, where <I> £ = cof(cr £ ), £/ie cofactor matrix of <J e :— D 2 u £ . Thus, using the 
divergence free property of the cofactor matrix <I> 6 (cf. Lemma 3.1) we can define the 
following alternative variational formulation for ( |2.1[ )-(2.2 1 : 



(<J £ ,fi)+ (div(n), Du £ ) = (g,u) 
(div(cr £ ),Dv) - ^-(<S> £ Du £ ,Dv) = (/ £ ,tj) 



en 



V/i G VF , 
Vtj g V . 



However, we shall not use the above weak formulation in this paper although it is 
interesting to compare mixed finite element methods based on the above two different 
but equivalent weak formulations. 



To discretize (2.5 1— ( 2.6 ) , let be a quasiuniform triangular or rectangular par- 
tition of if n — 2 and be a quasiuniform tetrahedral or 3-D rectangular mesh if 



() 
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n = 3. Let V h c be the Lagrange finite element space consisting of continuous 

piecewise polynomials of degree k(> 2) associated with the mesh T^. Let 



V g h := v h n v g , 

w h ,_ ^yh^nxn 



v> 1 :=v h nv , 

Wq := [V h ] nxn 



nWn 



In the 2-D case, the above choices of V and Wq are known as the Hermann- 
Miyoshi mixed finite element for the biharmonic equation (cf. j6j[16]). They form a 
stable pair which satisfies the inf-sup condition. We like to note that it is easy to check 
that the Hcrmann-Miyoshi mixed finite element also satisfies the inf-sup condition in 
3-D. See Section l3~2l for the details. 



Based on the weak formulation ( 2.5 1- ( 2.6 1 and using the above finite clement 



spaces we now define our Hcrmann-Miyoshi type mixed finite clement method for 



(Ol-OJ as follows: Find (u E h ,a%) E V g h x such that 



(2.7) 
(2.8) 



{o-%,Hh) + (div(nh),Du e h ) = (g,Hh) 
(dW(a e h ),Dv h ) + -(detK),^) = (f E ,v h ) 



v^h e ws 

Vv h E V h 



Let (a E ,u e ) be the solution to (2.5 (-(2.6 1 and (<r 



'in u h) 



solves ( 2.7 )-( 2.8 1. As 



u e is proved to be 
preserves the convexity even for 



mentioned in Section [T] the primary goal of this paper is derive error estimates for 
if — u e h and a e — a e h . To the end, we first need to prove existence and uniqueness 
of (<t|,u|). It turns out both tasks are not easy to accomplish due to the strong 
nonlincarity in (2.8 1. Unlike in the continuous PDE case where 
convex for all e (cf. |19]). it is far from clear if u e h 
small e and h. Without a guarantee of convexity for u £ h , we could not establish any 
stability result for u e h . This in turn makes proving existence and uniqueness a difficult 
and delicate task. In addition, again due to the strong nonlincarity, the standard 
perturbation technique for deriving error estimate for numerical approximations of 
mildly nonlinear problems does not work here. To overcome the difficulty, our idea 
is to adopt a combined fixed point and linearization technique which was used by 
the authors in [21) , where a nonlinear singular second order problem known as the 
inverse mean curvature flow was studied. We note that this combined fixed point and 
linearization technique kills three birds by one stone, that is, it simultaneously proves 
existence and uniqueness for ut and also yields the desired error estimates. In the 
next two sections, we shall give the detailed account about the technique and realize 



it for problem (2.7 1- (2 



3. Linearized problem and its finite element approximations. To build 
the necessary technical tools, in this section we shall derive and present a detailed 



study of the linearization of (2.5)-(2.6) and its mixed finite element approximations, 



First, we recall the following divergence-free row property for the cofactor matrices, 
which will be frequently used in later sections. We refer to [15, p. 440] for a short 
proof of the lemma. 

Lemma 3.1. Given a vector-valued function v = (ui,i>2) - '' ,"n) : f2 —> K™- 
Assume v E [C 2 (fl)] n . Then the cofactor matrix cof(-Dv) of the gradient matrix Dv 
of 'v satisfies the following row divergence- free property: 



(3.1) 



div(cof(Dv))i 



^Ta^cofpv)), 



fori = 1,2, 
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where (cof(£)v))j and (cof(_Dv))y denote respectively the ith row and the (i, j)-entry 
ofcof(Dv). 

3.1. Derivation of linearized problem. We note that for a given function w 
there holds 



&et(D z u E + tw) = det(L>V) + ttv(<$> e D 2 w) 



t n det(D 2 w). 



Thus, setting t — after differentiating with respect to t we find the linearization 



of M e (u e ) := -eA 2 u e + det(D 2 u E ) at the solution u e to be 



L u .(w) := -eA 2 w + tr($ £ D 2 w) = -sA 2 w + $ e 



D 2 w = -eA 2 w 



div($ £ Dw), 



where we have used (3.1 ) with v = Du £ . 

We now consider the following linear problem: 

(3.2) L u t(w) = q in Cl, 

(3.3) w = on 0Q, 

(3.4) D 2 w • v = on90. 



To introduce a mixed formulation for (3.2 1-( 3.4 1 , we rewrite the PDE as 

(3.5) X ~ D 2 w = 0, 

(3.6) -eAtr(x) + div($ e Dw) = q. 

Its variational formulation is then defined as: Given q G V *, find (x,w) G Wo x Vo 
such that 



(3.7) 
(3.8) 



(x,») + (dW(fi),Dw) = 

1 1 

(div(x),-Du) - -($ e Dw,Dv) 



V/i G W , 
V« G V . 



It is not hard to show that if (x, w) solves p?7]-(|3T8| then w G ff 2 (fi) n #o(ft) 
should be a weak solution to problem ( 3.2 1-( 3.4 1 . On the other hand, by the elliptic 
theory for linear PDEs (cf . [26 ) , we know that if q G V * , then the solution to problem 
(|3~2])-(|3~4"| satisfies w G H 3 (tt), so that x = D 2 w £ H l (Sl). It is easy to verify that 



(w,x) is a solution to (3.7)-([3 

3.2. Mixed finite element approximations of the linearized problem. 

Our finite clement method for ( 3.7 )-( 3.8 1 is defined by seeking (xh> w h) € x V 
such that 

(3.9) (Xfc.wO + (tiv(jM h ),Dw h ) = Wn h G Wtf, 

(3.10) (div( Xh ),Dv h ) - -(<S> e Dw h ,Dv h ) = (g,« h ) Vv h G itf. 



The objectives of this subsection are to first prove existence and uniqueness for 
problem (3.9 1-( 3.10 1 and then derive error estimates in various norms. First, we prove 
the following inf-sup condition for the mixed finite element pair (Wq , V Q h ). 

Lemma 3.2. For every Vh £ V h , there exists a constant (3q > 0, independent of 
h, such that 



(3.11) 



(dw(n h ),Dv h ) 
SU P ii — ii >Po|K||i?i. 
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Proof. Given v h € , set fi h = I nxn v h - Then (div(/.t /l ), Dv h ) = \\Dv h \\ 



L- 



> 



v h\\%i — PoWviiWh 1 Wf-hWii 1 ■ Here we have used Poincare inequality. □ 
Remark 3.1. By fTb\ Proposition 1], (3.11) implies that there exists a linear 
operator : W — > W h such that 

(3.12) (div(/i - U h p) , Dv h ) = Vv h e V h , 
and for fi G W n [iF(f2)] nXT \ r > 1, tfiere ZioZds 

(3.13) - n^Hijj < C/i'- J '||/x|| ff ! J = 0,1, l<I<min{* + l,r}. 

FKe note that the above results were proved in the 2-D case in }16\/ . However, they also 



hold in the 3-D case as (3.11 ) holds in 3-D. 

Theorem 3.1. For any q € Vq , there exists a unique solution (xh, Wh) € WqxVq 
to problem fiiM-Kwl. 



Proof. Since we are in the finite dimensional case and the problem is linear, it 
suffices to show uniqueness. Thus, suppose (xh,Wh) & Wq x Vq solves 



(Xh, fih) + (dW(n h ),Dw h ) = 
(div( Xh ),Dv h ) - -(&Dw h , Dv h ) = 

£ 



Vv h e V h . 



Let Hh — Xh, Vh — Wh and subtract two equations to obtain 

(Xh, Xh) + l -{<5> £ Dw hl Dw h ) = 0. 

Since u e is strictly convex, then $ e is positive definite. Thus, there exists > such 
that 

\\Xh\\h + - £ \\Dw h \\ 2 L2 <0. 

Hence, \h = 0, wt — 0, and the desired result follows. □ 

Theorem 3.2. Let (x,w) e [H r (ft)] nXn n W x ff r (») n V be the solution to 



(|377|)-(l3l3| and {xh,w h ) € W^ 1 x V^ 1 sofoes ([379]) -( [3T0| . TAen tfiere hold 

(3.14) ||x-Xh|U» ^Ce-t/i'-^HxIUi + lkllH'] 

(3-15) \\X-Xh\\m <Ce-ih l - 3 [\\x\\ H ' + \\w\\ H i] 

(3.16) Ik-tBhHffi ^Ce" 3 ^" 1 [\\x\\w + \\w\\ Hl ], 
Moreover, for k > 3 there also holds 

(3.17) \\w - w h \\ L 2 < Ce- 5 h l [\\ x \\ H ' + \\w\\h>]- 



Proof. Let IhW denote the standard finite element interplant of w in V$ . Then 

(3.18) (n ft x - Xh,^h) + (div(^ h ),D(I h w - w h )) 

= (IlfcX - X, Hh) + (div(nh),D(I h w - w)), 

(3.19) (div{IL h x-Xh),Dv h )- -(<S> £ D(I h w-w h ),Dv h ) 

£ 

= -~(<f> e D(I h w-w),Dv h ). 

£ 
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Let Hh = 11?,, — Xh and Vh = IhW — Wh and subtract (3.19) from (3.18) to get 

(n^x - Xh, n^X - Xh) + ~ ($ e D(I h w - w h ),D(I h w - w h )) 

= O^hX-X^hX-Xh) + (div(n ft x - Xh),D(I h w ~ w)) 



+ -(<f> £ D(I h w - w),D(I h w - w h )). 

e 



Thus, 



\\H-hX - Xhllia + - u>h)\\h 

< \\KhX- xlU 2 ll n /iX- X/iIIl 2 + l|n /l x-X/ l ||H 1 ll-D(^w-w) 



L 2 



+ -=-||i?(J fc «7 - to)|| ia ||2?(J h «; - «7 h )|| 



L 2 



< \\KhX - xh'W^hx - xhh* + Ch-^UhX - xhUAWhw ~ w)\\» 



+ -%\\D{I h w - w)\\ L 2\\D{I h w - w h )\\ 



L' 2 i 



where we have used the inverse inequality. 

Using the Schwarz inequality and rearranging terms yield 



\\RhX ~ Xh\\h + -\\D(I h w - w h )\\h 



(3.20) 

< C(||n h x - x\\h + h- 2 \\I h w - wf H1 + e- 3 \\I h w - wf H1 ). 
Hence, by the standard interpolation results [5) [TO] we have 

\\tthX-Xhh 2 < Cdln^x-xlU 2 + h- 1 \\I h w-w\\ H i +e-%\\I h w-w\\ JI i) 
<Ce-" 2 h l - 2 (\\x\\w + \H\ w ), 
which and the triangle inequality yield 

Hx-x.IIl 2 <Ce-ih l - 2 {\\ x \\ Hl + \\w\\ H ,). 

The above estimate and the inverse inequality yield 

llx - XftJIff 1 < llx - Ilfcxllff 1 + \\HhX - X/iIIh 1 

< llx - n h x||Hi + h-'WUhX - Xhh 2 
<Ce-ih l - 3 (\\x\\w + \H\w)- 



Next, from (3.201 we have 

\\D(I h w - w h )\\ L 2 < ^C[||n h x - xIIl 2 + h-'WDihw - w)\\ L 2 + e -^\\J h w - w\\ m ] 

(3.21) < Ce- 1 ^- 2 (||x||^ + II^IUO- 

To derive (3.161, we consider the following auxiliary problem: Find z £ H 2 {yi) n 
Hl(n) such that 

-eA 2 z + div(<[> £ Z)z) = -A (to - w h ) in CI, 
D 2 zv ■ v = on dn. 
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By the elliptic theory for linear PDEs (cf. ,26]), we know that the above problem has 
a unique solution z £ i?o(^) n H 3 (tt) and 

(3.22) \\z\\ H3 < C b (e)\\D(w-w h )\\ L 2 where C b (e) = O^ 1 ). 

Setting k — D 2 z, it is easy to verify that (k, z) £ Wo x Vq and 



+ (div(/x),Dz) = 



V/x e Wo, 



(div(«),Di;) - ^(<S> £ Dz,Dv) = ~{D(w - w h ),Dv) V« e V . 



It is easy to check that ( |3.9[ )-(3.10 1 produce the following error equations: 

(3.23) ( X -Xh^h) + (dW(fi h ),D(w-w h )) = V^G< 

(3.24) {dw{ X -Xh),Dv h )~- £ {^D{w~w h ),Dv h ) = {) W ft € V$. 



Thus, 



-11-0(10-™,,) Ilia = (dW(n),D(w-w h )) - l -{>$> £ Dz, D{w - w h )) 
div(K-U h K),D(w-w h )) - -(<f> e Dz,D(w - w h j) 



+ (&Vf{U h K),D{ W - w h )) 

div(« — U h K),D(w — Ihw)) — - ($ £ Dz, D(w — w h 

+ (xh - x,n fe K) 



)) 



div( K - U hK ), D(w - I h w)) - - ($ £ Dz, D(w - w h )) 

+ (xh - x, n^K - k) + (xh - x, k) 

div( K - U hK ), D(w ~ I h w)) - - (<S> £ Dz, D(w - w h )) 

+ (Xh - X, - k) + (div(x - Xh),Dz) 
div(K - Tl h n),D(w - I h w)) + (xh ~ X,^hK - k) 

+ (dW( x -Xh),D(z-I h z)) - -(<f> s D(w-w h ),D{z-I h z)) 

< \\div(K - n/,/c)||ia||£)(«J - hw)\\ L 2 + \\xh - x||i2||n fc « - k\\ L 2 
+ \\dW( X -Xh)\\L2\\D(z- hz)\\ L 2 



+ %\\D{z - I h z)\\ L 2\\D{w - w h )\\ 



L 2 



< c 



\\D(w - hw)\\ L 2 + h\\xh - xlU 2 + ^ 2 ||div(x - Xh)\\& 



+ \\D(w - w h )\\ L 2 \\z\\ 



H 3 ■ 



Then, by ( |3.14[ ),( |3.15| ),( |3.21| ), and ( |3.22| , we have 

\\D(w - w h )\\ L 2 < C b {e)e- 2 h l - l [\\ x \\w + \\w\ 



Substituting C b (e) = 0(e _1 ) we get (3.16). 
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To derive the L 2 -norm estimate for w — Wh, we consider the following auxiliary 
problem: Find (k,z) £ Wo x Vq such that 



(K,fi) + (div(n), Dz) = 



(div(«),Du) - - ($ e Dz, Dv) = -(to - to/j,u) 



V/i G W , 
Vv G V . 



Assume the above problem is 77 4 -regular, that is, z E H 4 (il) and 
(3.25) \\z\\ Hi < C b (s)\\w - w h \\ L 2 with C b {e) = 0{e- 1 ). 

We then have 

-\\w-w h \\ 2 L * = (div(K),D(w-w h )) ~ ^(<^ e D(w-w h ),Dz) 

= (div{U h K),D(w - w h )) - i ($ £ D(w - £>z) 
+ (div(« - 11/,/c), D(i£> - Wft) 

+ (div(« - U h K),D(w - hw)) 

= (xh - x, k) + (xh - x, n h K - «) 

- -($ £ £z, D(to - u> fc )) + (div(« - n hK ), £>(to - I fc u>)) 



1 



= (div(x - Xh), -D*) - -(* e ^(w - ^h), Dz) 

£ 

+ (Xh - X, n h K - k) + (div(« - n ft K), D(to - lira)) 
= (div( X -Xh),D(z-I h z))--(<f> s D(w-w h ),D(z-I h z)) 

£ 

+ (Xh - X, - k) + (div(« - Ilfe/c), D(to - J/jto)) 



< 



L 2 



[||div(x - Xh)IU» + ^\\D(w - w h )\\ L 2] \\D(z - I h z)\\ 
+ \\Xh - x\\hA\^hH - k\\ L 2 + ||div(«; - IL h n)\\ L 2 \\D(w - hw)\\ 

< Ch 3 [\\x~ XhWh 1 + - w h \\ H i]\\z\\ H i 

+ Ch 2 \\xh - xlU 3 ||«IU= + Ch\\w - I h w\\ H i\\K\\ H 2 
<C£- 5 h l (\\x\\w + \H\ w )\\z\\ Hi 

< CC b (s)E- 5 h l (\\x\\ H > + \\w\\ H i)\\w ~ W h \\ L 2, 



L 2 



where we have used (3. 14), (3. 151, (3. 16), (3.251, and the assumption k > 3. Dividing 
the above inequality by \\w — whWl 2 an d substituting C&(e) = 0(e~ 1 ) we get (3.171. 
The proof is complete. □ 



4. Error analysis for finite element method (2.7l-(2.8). The goal of this 
section is to derive error estimates for the finite clement method (2.7)-(2.8 1. Our main 
idea is to use a combined fixed point and linearization technique which was used by 
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the authors in [21] . 

Definition 4.1. Let T : x V g h -> x 6e a /wear mapping such that 
for any {p h ,v h ) £ W £ h x V g h , T{p h ,v h ) = (rW^, w h ), T( 2 )(^, w h )) safe/Zes 

(4.1) {fJ>h-T( lS> (fj,h,Vh),Kh) 

= (vth,Kh) + (drv(Kh),Dvh) - (g,Kh) V/c/, e Wq , 

(4.2) (div^-TW^,^)),^) - i^^K-T^^,^)),^) 

= (div(/i ;i ),Dz ft ) + ~(det((j,h),Zh) ~ (f £ i z h) Vzh G Vq. 



By Theorem 3.1 we conclude that T(ph,Vh) is well defined. Clearly, any fixed 
point (xh,Wh) of the mapping T (i.e., T(xh,Wh) — (Xh,Wh)) is a solution to problem 
(2.7)-(|2.8|), and vice- versa. The rest of this section shows that indeed the mapping 



T has a unique fixed point in a small neighborhood of (Ih& e , Ih uE )- To this end, we 
define 

B h {p) := {{p h ,v h ) € x V g h ; \\p h - I h a £ \\ L 2 + -±=\\v h - I h u £ \\ m < p}. 

v £ 

Z h := {(p h ,v h ) e x V/ 1 ; (/i h) K h ) + (div( Kh ), Dv h ) = (g,K h ) V Kfl £ W*}. 
B h (p) :=B h (p)nZ h . 

We also assume a e S H r (il) and set Z = min{fc + 1, r}. 

The next lemma measures the distance between the center of Bh(p) and its image 
under the mapping T. 

Lemma 4.1. The mapping T satisfies the following estimates: 

(4.3) \\I h a £ -T^(I h a £ ,I h u £ )\\ m < C 1 (e)h l - a [\\o'\\ H > + \\u e \\ H i], 

(4.4) \\I h a £ -T^(I h a £ ,I h u £ )\\ L2 < C 2 (e)h l ~ 2 [\\a s \\ Hl + ||u s || ff! ] , 

(4.5) \\hu s -T^ (I h a s ,I h u s )\\ H i < C^h^lW^WH, + \\u s \\ Hl ]. 



Proof. We divide the proof into four steps. 

Step 1: To ease notation we set u>h — Ih<3 e — T^(Iha e ,IhU e ), Sh = IhU e — 
T^{I h a £ ,I h u e ). By the definition of T we have for any [ph, Vh) € Wq x V h 

(uJh,Ph) + (div(p h ),Ds h ) = (I h (j e ,p h ) + (div(p h ),D(I h u s )) - (g,Ph), 
(div{w h ),Dv h ) - ^(<S> e Ds h ,Dv h ) = (div(I h a £ ),Dv h ) + l -(det{I h a e ),v h ) - {f,v h ). 



It follows from p^J-p^J that for any (ph,v h ) € Wfi x V h 



(4.6) ((dh,lih) + (div(p h ),Ds h ) = (I h a E - a e ,p h ) + (div{p h ), D(I h u E - u e )), 

(4.7) (div(w h ), !>»/») - ^(® s Ds h ,Dv h ) = (div(I h a £ -a £ ),Dv h ) 

+ -(det(I h a e )-det(a £ ),v h ). 
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Letting vh — s/i, fJ-h = m (4.6)-(4.7), subtracting the two equations and using 
the Mean Value Theorem we get 

(u h ,u h ) + -^ £ Ds h ,Ds h ) = (I h a £ -a e ,uj h ) + (div(w fc ), D(I h u e - u £ )) 

+ (div(<r - I h <r E ),Ds h ) + -(det(<r £ ) - det(I h a e ),s h ) 
= (I h a e - a e ,u h ) + (div(u h ),D(I h u e - u e )) 
+ (div(a - I h a s ),Ds h ) + j(* e : (a £ - I h a E ), s h ), 

where * e = coi(t//j<7 e + [1 - r](T £ ) for r G [0, 1]. 

5^ep 2: T/ie case n = 2. Since $ E is a 2 x 2 matrix whose entries are same as 
those of rIh<J E + [1 — t]ct £ , then by ( |1.11[ ) we have 

\\V S \\ L 2 = ||coi(t/^ £ + [1 - T]a e )\\ L 2 = \\ T I h a £ + [1 - T ]a s \\ L 2 
< \\I h a e \\ L s + \\<t s \\ L 2 < C\\a e \\ L 2 = 0(s~i). 

Step 3: The case n = 3. Note that (* E )ij = (cof(Tl h cr £ + [1 - r]a £ ))ij = 
det(r/^(T e |j : j + [1 — rjo-^ljj), where cr e \ij denotes the 2x2 matrix after deleting the 
ith row and jth column of a e . We can thus conclude that 

|(* e ) y | < 2 max (\r(I h a s ) st + [1 - r](<7 e ) flt |) 2 
<C m^.\(*% t \ 2 <C\\a e \\ 2 Loa . 

Thus, ( |1.11D implies that 

||* e ||L* <C\\a s \\l aa =0(e- 2 ). 
Step 4- Using the estimates of ||^ £ ||l 2 we have 

Q 

\\u h f L 2 + -\\Ds h \\ 2 L2 < \\I h a £ - a e \\ L 2\\uj h \\ L 2 + ||«h||Hi||I>(Zj,u e - u e )\\ L 2 

+ \\hv E - cr e \\ m \\Ds h \\ L 2 + C(e)\\a E - I h a e \\ m \\s h \\ m , 

where we have used Sobolev inequality. It follows from Poincare inequality, Schwarz 
inequality, and the inverse inequality that 

(4.8) \\w h \\ 2 L . + 6 -\\s h \\ 2 m < C{e)\\I h <r e - a e \\ 2 H1 + C\\uj h \\ m \\I h u £ - u s \\ m 

< C(s)h 2l - 2 \\^\\ 2 Hl + Ch- l \\u h \\ L 4I h u E - u s \\ H ,. 



Hence, 



Therefore, 



+ ~KHk < C(e)h 2l - 2 \\a*\\ 2 Hl + Ch 2l -*\\u E \\ 2 Hl . 



\\Lo h \\ L 2<C 2 {e)h l - 2 [\\a E \\ Hl + \\u E \\ Hl ], 
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which and the inverse inequality yield 

IKIk^Ci^ft'-^iKH^ + H^iiH,]. 

Next, from ( |4.6[ ) we have 

(div(fi h ),Ds h ) < \\uj h \\ L 2\\fi h \\ L 2 + \\I h a e - a e \\ L 2 \\fi h \\ L 2 
+ \\div(^ h )\\ L 2\\D(I h u E - u £ )\\ L 2 



<C 2 (s)h l - 2 [\\a s \\ Hl + \\u s \\ Hl ]\\fi h \\ 



H 1 - 



It follows from (3.11 1 that 

\\Ds h \\ L ,<C{e)h l - 2 [\\o*\\ H i + \\u*\\ H i] 



(4.9) 



To prove (4.5), let (k,z) be the solution to 
(«,/*) + (div((i),Dz) = 



(div(K),Du) - -(<S> e Dz,Dv) = -(Ds h ,Dv) 

e e 



Vm g Wo, 

Vv G Vb, 



|*||*» < c 6 (£)||-D S/l || 



L 2 • 



and satisfy 



Then, 

^\\D Sh \\ 2 L 2 = (div(/c),Ua h ) - ~($ e Dz,Ds h ) 

= (div(n^),^ S/l ) - ~($ e x>*,£>« h ) 

= -(w h ,n h «) - i($ B D«, Ds h ) + (I h a £ - a e ,U h n) 



+ (div(U h n),D(I h u s -u £ )) 
= -(oj h , k) + (oj h , k - UhKj - -($ e Dz, Ds h ) 

£ 

+ {I h a e - a £ ,U h K) + (div(U h K),D(I h u E - u e )) 

= (div((j h ),Dz) - -($ e Ds h , Dz) + (uj h , k - IL h K) 
e 

+ (I h a £ - <r E ,Il hK ) + (div(U h K),D(I h u s - u e )) 
= (div(uj h ),D(z ~ I h z)) ~ -($ e Ds h , D(z - I h z)) + (u h) k - U h n) 

£ 

+ {I h a £ - a £ ,\l h n) + (div(Tl h K),D(I h u e - u 6 )) 

+ (div(a £ - I h a £ ),I h z) + -(det(a £ ) - det(I h a e ), I h z) 

£ 

< ||div(w h )IM|£(* - I h z)\\ L 2 + ~\\<Z> £ \\ La o\\Ds h \\ L 2\\D{z - I h z)\\ L 2 

+ II^/iIU 2 II k - n^n^ + \\i h cr £ - cr e || L 2||n, l K|| L 2 

+ ||div(n h K)|| L »||D(/ h u e -« 6 )|| La 

+ ||div(a e - I h a e )\\ L 2\\I hZ \\ L 2 + -H^ll^ll^ - I h a e \\ m \\I h z\\ m 

£ 

<Ch 2 (\\ u\\m + -^\\Ds h \\ L 2)\\z\\ H 3 + C(£)h l 1 (\\I h z\\ L 2 + \\I h z\\ H i)\\c 



Ch\\u h \\ L 2\\ K \\ m +Ch l \\a e \\ H ,\\U hK \\ L 2+Ch l - i \\U h K\\ H i\\u s \\ H , 
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<C 2 ( e )e- a fc , - 1 [||u e || ff , + ||(7 e || ff .]||z||H3 
<C 2 (e)e- a C 6 (e)/i , - 1 [||u e || H , + ||(7 e || H ,]||D« h || ia . 
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Dividing by H-Ds/Jl^a, we get (4.51. The proof is complete. □ 

Remark 4.1. Tracing the dependence of all constants on e, we find that C\{e) = 
0(1), C 2 (e) = 0(1), C 3 (e) = 0(£~ 2 ) when n = 2, and d(e) = 0{eT^), C a (e) = 
0(e~5) ; (73(e) = 0(e~5) wften n = 3. 

The next lemma shows the contractiveness of the mapping T. 

Lemma 4.2. There exists an ho = o(e^) and po — o(e 13 1 log h\ n ~ 3 h%~ 1 ), such 
that for h < ho, T is a contracting mapping in the ball Bh(po) with a contraction 
factor |. That is, for any (/ift,Uft), (XhiWh) £ Bf t (p ) there holds 



(4.10) iitW^^^-tW^^iu. + ^iit^^,^)-^ 2 )^,^)!! 



H 1 



Proof. We divide the proof into five steps. 
5iep i: To ease notation, let 

TV =TW(p Jh ,v h )-TW(xh,w h ), T^ =T^(p h ,v h )-T^( Xh ,w h ). 
By the definition of T^ % > we get 

(4.11) (T^,K h ) + (dW( Kh ),D(T^)) = V/s h eW h , 

(4.12) (div(T«),£>* h ) - l($ e i?(rW),Dz h ) 

= - [($ 8 Z>(u; A - v h ),Dz h ) + (det( Xh ) - det(p h ),z h )] Vz h e V h . 

Letting z/, = T^ 2 ' and «/, = T' 1 ), subtracting ( |4.12 ) from (4.11 ), and using the Mean 
Value Theorem we have 

(T^\T {l) ) + -(^ E DT {2) ,DT {2) ) 

£ 

= l - [(<f> E D(v h - w h ),DT&>) + (det(p h ) - det(xh), T (2) )] 
- - [($ e Z?(« fc - w h ), DT&>) + (A h : (p h - T«)] 
= 1 [(&D(v h - w h ),DT&>) + (* £ : (p h - Xh ), T (2) ) 

£ 

+ ((A h ~^):(p h - Xh ),T^) 
= ^[(div(^T^),D(v h -w h )) + 

+ ((A h -&):(p h - Xh ),T^) 
= -[(div(n fe ($ e T( 2 ))), D(v h -w h 



) + (p h - Xh ,^T^) 



+ ((A h -^):(p h - Xh ),T^)} 

1 [(d> e r( 2 ) - U h (^T^),p h - Xh) + ((Ah - <& E ) : (w, ~ Xh), T (2) )] 
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< 1[||<FT( 2 ) -U h (^T^)\\ L 4f, h - xhh* 

+ C\\A h -^\\ L 2\\^ h - Xh \\ L2 \\T^\\ L ^] 

< j[\\9 e T^-U h ^ e T^)\\ L 4^ h - X h\\L' 

+ \iogh\ 3 - n h 1 - 2 i\\A h -^\\ L 4^ h -x h \\L4T {2) \\H^, 

where = cof(fih + T~(xh — M/i))> T € [0, 1]. n = 2,3. We have used the inverse 
inequality to get the last inequality above. 

Step 2: The case of n — 2. We bound ||<I> e - A h \\ L 2 as follows: 

||$ £ - K h \\ L 2 = ||cof(a £ ) - cof{fi h + T( Xh - Hh))\\L* 
= \W £ ~ Vh~ r(xh - A*h)IU 2 

< ||cr e - I h (7 £ \\ L 2 + \\I h <7 £ - HhWv + \\Xh- HhWl? 

<Ch l \\a £ \\ Hl +3p . 

Step 3: The case of n = 3. To bound ||<I> e — A^Wl 2 in this case, we first write 

||($ £ - A h )ij\\ L 2 = ||(cof(o- e )ij) - co{{/i h + r(xh - 

= ||det(cr e |^) -det(ii h \ij + i~(Xh\ij - Wi|ij))IU 2 > 

where a\ij denotes the 2x2 matrix after deleting the i th row and j th column. Then, 
use the Mean Value theorem to get 

||($ £ - K h )ij\\ L 2 = ||detO £ | 4J ) - det(fj, h \ij + r(xh\ij - Mfc|ij))IU 2 
= ll A y ■ (o- e \ij - Hh\ij ~ r(xh\ij ~ M/i|ij))IU 2 
< ||Ajj||z,oo \\cr E \ij - Hh\ij - r(xh\ij - Wi|i?)IU 2 ; 

where Ay = cof(<7 E |ij + X(fJ,\ij - r(xh\ij - lAij) - v £ \ij)), X € [0, 1]. 
On noting that Ay e R 2 , we have 

IIAij-Hioo = ||cof((T £ |y + X{ii\ij - r{xh\i 3 - lAa) - v e \ij))\\L°° 

= \W E \%3 + Hl^Uj - T (Xh\ij - - 0- E \ij)\\ L co 

<CMl«<j. 

Combining the above estimates gives 

C 

||($ £ - A h )ij\\^ < -h%j - Mhlij - r(xh\ij - fJ-h\ij)\\Li 

<j(h l \\^\\ Hl +p ). 
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Step 4: We now bound ||$ £ T( 2 ) - n /l ($ £ T^)|| L 2 as follows: 



H 1 



= C/i 2 ( 

< C7i 2 ( 

< Ch 2 ( 

< Ch 2 ( 

< Ch 2 ( 
Ch 2 



< 



13 
£ 6 



&TW\\l 2 + \\D(<f> e TW)\\l 2 ) 

^T^\\\ 2 + \\^ £ DT^\\ 2 L2 + \\D<5> E T^\\ 2 L2 ) 

$ e ||i=o||DT (2) |||a 1 lln ^ £|12 



$ e || 2 4 ||T^| 



2 

L 1 



D$ s ||i3||T( 2 )||| 6 ) 



<$> e \\ 2 tA\tW\\ 2 h1 + ll^lllocljOT^H 2 , + ||M e ||| 3 \\T^f m ) 



L 4 I 



|D$ £ || 2 3 )||Dr( 2 )|| 



2 

L 2 



\\DTW\\ 



L 2 ' 



where we have used Sobolev's inequality followed by Poincare's inequality. Thus, 

Ch, 



L 2 — lis 

£12 



||i?T^|| 



L 2 ■ 



Step 5: Finishing up. Substituting all estimates from Steps 2-4 into Step 1, and 
using the fact that <I> e is positive definite we obtain for n — 2, 3 



||T«|| 2 2 + -IIOT^H 2 , < Ce-v(h+ llogh^h 1 -^ po)\\n h - xh\\w 



\\DTW\ 



L 2 ■ 



Using Schwarz's inequality we get 

l|T (1) || L 2 + ^||T( 2 )|| ff i < Ce-$(h+ llog/il 3 ""/! 1 -" 



POlWVh - Xh\\L 2 - 



Choosing ha = o{e" 2 ) and po — o(e^ | log/i| n 3 h 2 1 ), then for h < ho there 
holds 

||t«iu 3 + ^||t( 2 )|| h1 <~y h -xh\\» 

< ^(llMft - XftlU 2 + ^W Vh ~ w h\\m)- 

The proof is complete. □ 

We are now ready to state and prove the main theorem of this paper. 

Theorem 4.1. Let Pl = 2[C 2 {£)h 1 - 2 + g ^-h l - 1 ]{\\a £ \\ H i + \\u e \\ H i). Then there 
exists an hi > such that for h < min{hQ, hi}, there exists a unique solution {o~ e l ,u e h ) 
to (|2~7)-(|2~I| in the ball B h (pi). Moreover, 

(4.13) - + -^-\\u' - u%\\ m < C 4 (e)h l - 2 (\\<j'\\ Hl + \\u e \\ Hl ), 



(4.14) 



W e - v e h \\m < C 6 (e)/i , - 3 (||<7 b ||hi + KM 



Proof. Let Qu/i,i>/i) € Bh(pi) and choose hi > such that 

2 

/ 25 \ 2!-re 

2(3-n) / £12 \ 

/ii| log /ii| 21 -™ < C - . ... — n n — n — r and 

" \C' 3 (e)(\\^\\ Hl + \\^\\ Hl ) 



2(3-n) / £ 12 

/ii| log /ill 2 '""" 2 < C 



C 2 (e)Q\ae\\ H i + \\ut\\ Hl ) 
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Then h < minj/io, hi} implies p\ < pq. Thus, using the triangle inequality and 
Lemmas |4. 1 1 and |4. 2| we get 



+ \\TW(I h a s ,I h u e )-TM(» h ,v h )\\ L2 + -^\\hu s -T^(I h a £ , I h u s )\\ H i 

V £ 

+ -^\\TW(I h a e ,I h u s ) - TW(» h ,v h )\\ m 

<[c 2 ( £ )/i'- 2 + ^/ l '- 1 ](|K|| H , + K|| H o 

s P 1 , P 1 ^ 1 

< Y + T = P!<1. 

So T(/j,h,Vh) € Bh(pi). Clearly, T is a continuous mapping. Thus, T has a unique 



fixed point (cr^,u|) € Bh{pi) which is the unique solution to ( 2.7 1-( 2.8 1. 
Next, we use the triangle inequality to get 

\W e - + ^\W e - < h e - h* e \\» + \\h<r £ - o%\\ L * 

v £ 

+ "pOK ~ ^"Iff 1 + 11-^ - u V\\m) 
v £ 

<Pi + Ch l - 1 (\\a'\\ Hl + \\u e \\ H i) 
<C 4 (e)h l - 2 (\\a E \\ Hl + \\u s \\ Hl ). 
Finally, using the inverse inequality we have 

h £ - v E h \\m < K - h<T E \\m + \\h° e - <7 E h \\m 

< IK - h<T £ \\ H i + Ch^\\I h a £ - a e h \\ L 2 

< Ch l - l \\a e \\ H i + Oh' 1 pi 
<Cs{e)h l -%<j e \\ H i + \\u s \\ Hl ]. 

The proof is complete. □ 

Remark 4.2. By the definition of p\, and the remark following Lemma \4-1\ we 
see that C^{e) = Cs(e) = 0(e~%) when n = 2, C±{s) — Cs(e) = 0(e -4 ) when n = 3. 

Comparing with error estimates for the linearized problem in Theorem |3.2| we 
see that the above iJ 1 -error for the scalar variable is not optimal. Next, we shall 
employ a similar duality argument as used in the proof of Theorem |3.2| to show that 
the estimate can be improved to optimal order. 



Theorem 4.2. Under the same hypothesis of Theorem J^.l there holds 



(4.15) ||u e - u%\\ H i < C 4 (e)e- 2 [h 1 ' 1 + C 5 (e)h 2 ^] (|K|| H , + \\u e \\ H >). 

Proof. The regularity assumption implies that there exists («, z) £ Wo X Vq H 
H 3 (n) such that 

(4.16) (k, fjt) + (div(p), Dz) = Wp e W , 

(4.17) (div( k),Dv) - -(<S> e Dz,Dv) = -(D(u e - ul),Dv) Vw G V , 

e e 
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with 



(4.18) 



\z\\H><C b ( S )\\D(v*-u° h )\\ 



L 2 ■ 



It is easy to check that a £ — a\ and u £ — ul satisfy the following error equations: 

(4.19) (a £ - a e h , + (div(/i ft ), D{u £ - u%)) = V/x ft e Wtf , 

(4.20) (div(a £ - 4),£>0 + -(det(a e ) - det(^),« fc ) = Vu h € V^ 1 . 



By (4.16 1-(4.20 ) and the Mean Value Theorem we get 



l -\\D{u £ -u £ h )\\l 2 = (dxv(K),D(u s -ui))-±(&Dz,D(u s -u%)) 

= (div(U h K),D(u s -u%)) - ^($ £ D(u £ -u%),Dz) + (div(n~Il h K),D(u e - u £ h )) 
= (p% - <r £ , U h k) - J (<f> £ D(u £ ~ u £ h ),Dz) + (div(« - U hK ), D(u £ - u%)) 

= (a% - a £ , K ) - \ ($ £ D(u £ - u £ h ),Dz) 

+ (div(n -Tl h K),D{u £ -I h u £ )) + (a £ h - a £ ,Jl h K- k) 
= (div(<7 E - a%),Dz) - -^ £ D(u £ - u £ h ), Dz) 

+ (div(« - IL h K), D{u £ - I h u £ )) + (a £ h - a £ ,Tl h K - k) 
= (div(a s -a s h ),D(z-I h z)) - D(u £ - u £ h ),D{z ~ I h z)) 

+ (div(/c - IL h K), D(u £ - I h u £ )) + (a% - a £ ,U h n - k) 

- i(det(a £ ) - det(ag), I h z) - - £ {<S> £ D(u £ - u £ h ),D{I h z)) 

= (div(a £ - a%),D(z - I h z)) - -^ £ D{u £ - u%),D(z - I h z)) 
+ (div(« - IL h K), D(u £ - I h u £ )) + (a £ h - a £ ,Tl h K - k) 

- : (<x £ - ^),/ /l z) - ^($ e U(u e - u s h ),D(I h z)), 

where \& e = cof(cr e + t[g\ - ct £ ]) for r e [0, 1]. 
Next, we note that 



(# E : (a £ - a £ h ),I h z) + (<!> £ D(u £ - u e h ),D(I h z)) 

= ($ £ : (a £ - a £ h )J h z) + (div($ £ I h z), D(u £ - u £ h )) + ((* £ - $ £ ) : (a £ - a £ h )J h z) 
= (a £ - a s h ),$ s I h z) + (div(Il h (<Z> £ I h z)),D(u s - u £ h )) + ((* £ - $ £ ) : (a £ - a £ h ),I h z) 

+ (div($ £ 4z - Tl h ($ E I h z)), D(u £ - I h u £ )) 
= (<r £ - a £ h ^ £ I h z - II h (& a I h z)) + ((^ £ - $ £ ) : (a £ - <J £ h ),I h z) 

+ (div(<Z> e I h z - IL h (<Z> s I h z)),D(u e - I h u £ )). 
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Using this and the same technique used in Step 4 of Lemma |4.2| we have 

-JD(u e - u%)\\%, = (div(<7 £ - a%),D{z - I h z)) - - £ (<S> E D{u E - u E h ),D(z - I h z)) 

+ X - [(($ £ - : (a s - a e h ),I h z) + (a E - a s h ,Il h (<Z> s I h z) - <S> E I h z) 

+ (div{Il h ($ s I h z) - <5> E I h z),D{u E - 4u e ))] + (<j e h - a e ,U hK - K ) 
+ (div(« - H h n) , D(u E -I h u E )) 
< [||div(a £ - aDWv + ~\\D(u e - ul)\\ L2 ] \\D(z - I h z)\\„ 

+ _ ^|| L2 ||a e - afJI^H^zll^ + \\a E - a%\\^\\n h {^I h z) - $ e I h z\\ L 2 

+ \\div(Il h ($ e I h z) ~ $ e I h z)\\ L 4D(u s - I h u E )\\ L2 ] +\\k- Il h K\\ L 2\\a E - al\\ L 2 
+ ||div(« - Il hK )\\ L 2\\D(u E - I h u E )\\ L 2 



< Ch 2 (\\a E - a E h \\ m + -^\\u E - u E h \\ m )\\ Z \\ H3 



C 



< 



+ Ch\\a E - a E h \\ L 2\\ K \\ H i + C\\u E - I h u E \\ m \\k\\ H i 
f (Ci(e) + C 5 ( e))h 



— [IKIIw + ll«'IU<] + °' ( % h ' 2 ll«* - *'IU-}W"" 



<aM{ (c ' M+ .? ( " )ft '" [ii^h: + iKi| g .] 



£2 



C 4 (e)/i 



;-2 



* E IU*}lP( 



We now bound ||<I> e — \J/ e 1 1 ^2 separately for the cases n — 2 and n — 3. First, 
when n = 2 we have 

||<i> £ - || L2 = ||cof(<7 E ) - cof(a| + r[a E - a%])\\ L 2 

= h e -{*t + T[<T E -*l])\\ L * 

<C 4 (s)h l - 2 [\\a E \\ Hl + \\u E \\ Hl ]. 

Second, when n = 3, on noting that 

|($ e - *%-| = \{coi{a E )) l3 - (cof(4 + T[a £ - a E h ])) l3 \ 

= \det(a e \ij) - det(a E \ tJ + r[a E \ l3 - <^| tf ])|, 

and using the Mean Value Theorem and Sobolev inequality we get 

!!(*%■ " mnWv = (1 - t)||(A^ : (a E \ l3 - a E h \ l3 )\\ L 2 

<\wr\\ H A\^k-<v 3 \\m, 

where (A e )^' = cof(o- e |y + \[o%\ij - a E \ l3 ]) for A e [0, 1]. Since (A e )« € R 2x2 , then 



Thus, 
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Finally, combining the above estimates we obtain 

\\D(u £ - ul)\\ L , < C^e- 2 [h 1 - 1 + C^h 2 ^] (||a e || ff! + \\u*\\ H ,). 

We note that 2(1 - 2) > I - 1 for k > 2. The proof is complete. □ 

5. Numerical experiments and rates of convergence. In this section, we 
provide several 2-D numerical experiments to gauge the efficiency of the mixed finite 
element method developed in the previous sections. We numerically determine the 
"best" choice of the mesh size h in terms of e, and rates of convergence for both u° — u B 
and u £ — u\. All tests given below are done on domain — [0, l] 2 . We refer the reader 
to [19, 28 for more extensive 2-D and 3-D numerical simulations. We like to remark 
that the mixed finite element methods we tested are often 10-20 times faster than the 
Aygris finite element Galerkin method studied in [20] . 

Test 1:. For this test, we calculate — wfj| for fixed h — 0.015, while varying e 
in order to estimate \\u e —u \\. We use quadratic Lagrange element for both variables 



and solve problem (2.51(2.6) with the following test functions: 



1 2,2 2,2 2,2 

(a) . u" = -e i , f = (1 + x 2 + y*)e * , g = e = , 

(b) . u° = x 4 + y 2 , f = 2Ax 2 , g = x 4 + y 2 . 



After having computed the error, we divide it by various powers of e to estimate 
the rate at which each norm converges. Tables 5.2 and 5.4 clearly show that \\<r Q — 
"Ilk 2 — 0{ £Z )- Since h is very small, we then have - 



\H 2 



(9(e<i). Based on this heuristic argument, we predict that — u 6 \\ H 2 = 0(e*). 

0{e) and \\u°-u e \\ H i « 



Similarly, from Tables 5.2 
0( £ s). 



and 



5.4 



we see that ||u — it e || ^ 2 





FlG. 5.1. Test la. Computed solution u e h (left) and its L 2 -error (right) (s — 0.05) 



Test 2:. The purpose of this test is to calculate the rate of convergence of \\u e — 
for fixed e in various norms. We use quadratic Lagrange element for both variables 



and solve problem (2.5 1 — (2.6 > with boundary condition D u e v 



e on dVt being 
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£ 


ii =■ (in 

K-« U IU" 


IK -"Iff* 


II cr nil 


0.75 


0.031968735 


0.168237927 


1.412579201 


0.5 


0.038716921 


0.196397556 


1.559234748 


0.25 


0.040987803 


0.206004854 


1.644877503 


0.1 


0.032218007 


0.168139823 


1.541246898 


0.075 


0.028113177 


0.150389494 


1.480968264 


0.05 


0.022258985 


0.124863926 


1.386775396 


0.025 


0.013676045 


0.086203248 


1.217747100 


0.0125 


0.007816727 


0.057280014 


1.052222885 


0.005 


0.003511072 


0.032109189 


0.853140082 


0.0025 


0.001863935 


0.020252025 


0.722844382 


0.00125 


0.000973479 


0.012568349 


0.611218455 


0.0005 


0.000404799 


0.006544116 


0.492454059 



Table 5.1 



Test la: Change of \\u° - u% || w.r.t. s (h = 0.015J 



e 


IK-« u llta 


\\ul-u u \\ H i 


lkj-o- u || l2 


£ 






0.75 


0.04262498 


0.194264425 


1.517915136 


0.5 


0.077433843 


0.277748087 


1.854253057 


0.25 


0.163951212 


0.412009709 


2.326208073 


0.1 


0.322180074 


0.531704805 


2.740767624 


0.075 


0.374842355 


0.54914479 


2.829960907 


0.05 


0.445179694 


0.558408453 


2.932672906 


0.025 


0.54704179 


0.545197212 


3.062471825 


0.0125 


0.625338155 


0.51232802 


3.146880418 


0.005 


0.702214497 


0.454092502 


3.208321232 


0.0025 


0.745574141 


0.405040492 


3.232658349 


0.00125 


0.778783297 


0.355486596 


3.250640603 


0.0005 


0.809598913 


0.29266175 


3.293238774 



Table 5.2 



Test la: Change of \\u° - u% || w.r.t. e (h = 0.015) 



replaced by D 2 u e v ■ v — h £ on dil and using the following test functions: 



(a) . U z = 20x 6 + y 6 , f 

g £ = 20x 6 + y 6 , 

(b) . u £ = ssin(x) + ysia(y), f 



„4„,4 



e(7200a; 2 + 360?/ 2 ), 



h E = 600x 4 " 2 



30 y 4 ^ 2 . 



g e = a;sin(a;) + ysin(y), h £ 



lSOOOa; 4 ^ 

v. 

(2cos(ir) — xsin(a;))(2cos(i/) — y * sin(y)) 

— e(xsm(x) — 4cos(a;) + ysin(y) — 4cos(y)) 
(2cos(ie) — a;sin(x))i/ 2 + (2cos(y) — ysm{y))v. 



After having computed the error in different norms, we divided each value by a 
power of h expected to be the convergence rate by the analysis in the previous section. 
As seen from Tables 5.6 and 5.8 the error converges exactly as expected in ff 1 -norm, 
but erf appears to converge one order of h better than the analysis shows. In addition, 
the error seems to converge optimally in L 2 -norm although a theoretical proof of such 
a result has not yet been proved. 



MIXED FINITE ELEMENT METHODS FOR MONGE-AMPERE EQUATIONS 



23 



£ 


II cr fill 

\K-u u \\ L 2 


K-« u ||jn 


II c- fill 


0.75 


0.080523289 


0.441995475 


3.65931592 


0.5 


0.082589346 


0.448160685 


3.706413496 


0.25 


0.074746237 


0.412192916 


3.603993202 


0.1 


0.051429563 


0.309140745 


3.233364656 


0.075 


0.043554563 


0.273452007 


3.091264143 


0.05 


0.033436507 


0.226024335 


2.885806127 


0.025 


0.020115546 


0.158107558 


2.538905473 


0.0125 


0.011590349 


0.107777549 


2.211633785 


0.005 


0.005376049 


0.06303967 


1.820550192 


0.0025 


0.002939459 


0.041182521 


1.559730105 


0.00125 


0.001580308 


0.026467488 


1.330131572 


0.0005 


0.000679181 


0.014385878 


1.075465946 



Table 5.3 



Test lb: Change of \\u° - u%\\ w.r.t. e (h = 0.015J 





IK-« u lli,2 


i s On 




£ 


e 






0.75 


0.107364385 


0.510372413 


3.932190858 


0.5 


0.165178691 


0.63379492 


4.4076933 


0.25 


0.298984949 


0.824385832 


5.096816065 


0.1 


0.514295635 


0.977588871 


5.749825793 


0.075 


0.580727513 


0.99850555 


5.907052088 


0.05 


0.668730140 


1.010811555 


6.102736940 


0.025 


0.804621849 


0.999959999 


6.385009233 


0.0125 


0.927227955 


0.963991701 


6.614327771 


0.005 


1.075209747 


0.891515564 


6.846366682 


0.0025 


1.175783722 


0.823650411 


6.975325082 


0.00125 


1.264246558 


0.748613599 


7.074033284 


0.0005 


1.358362838 


0.643356045 


7.192074244 



Table 5.4 



Test lb: Change of \\u° - u||| w.r.t. e (h = 0.015,) 



Test 3. In this test, we fix a relation between e and h, and then determine the 
"best" choice for h in terms of e such that the global error u° — u £ h has the same 
convergence rate as that of u° — u £ . We solve problem (2.5 1— ( 2.6 ) with the following 
test functions: 



(a) . u° = x 4 + y 2 , 

(b) . u° = 20x 6 + y 6 , 



f = 18000x 4 y 4 , g = 20x 6 + y e . 



To see which relation gives the sought-after convergence rate, we compare the 
data with a function, y = (3x a , where a = 1 in the L 2 -case, a = \ in the -H^-case, 
and a = 4 in the H 2 -case. The constant, (3 is determined using a least squares fitting 
algorithm based on the data. 

As seen in the figures bel ow, the best h — e relation depends on which norm one 

0(e) 
e, 



considers. Figures 5. 3j|^4]and[5j][5!8] indicate that when h = e 2 , 
and ||<7° — ct^Wl 2 ~ 0(ei ). It can also be seen from Figures 5.5 - 5.6 that when h 



l h\\L 2 
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FlG. 5.2. Test la. Computed Solution and its L -error (h — 0.015J 



n 


h 


\\u £ -U%\\ L 2 


\W e -u £ h \\ H i 


lk e -°1IU 2 




10 


0.1 


0.004334849 


0.335913679 


0.083695878 


5.995796194 


20 


0.05 


0.000545090 


0.084457090 


0.011891926 


1.813405912 


30 


0.033333333 


0.000161694 


0.037576588 


0.003840822 


0.916912755 


40 


0.025 


6.82423E-05 


0.021145181 


0.001747951 


0.574128035 


50 


0.02 


3.49467E-05 


0.013535235 


0.000959941 


0.403471189 



Table 5.5 

Test 2a: Change of \\u c - u%\\ w.r.t. h (e = 0.001,) 



\\u -u%\\ m =O(si). 
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n 


h 


IK-^hlLa 


IK—"hllffi 


ll^" -0 * II £2 


h 3 


h 2 




10 


0.1 
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20 


0.05 
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30 
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Test 2a: Change of \\u c - u%\\ w.r.t. h (e - O.OOlj 



n 


h 


IK - "IIU 2 


~ u I\\hi 




\W e ~°l\\m 


10 


0.1 
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30 
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4.99964E-07 
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40 
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50 


0.02 


1.07997E-07 
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1.20594E-06 


0.000553558 
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ll«*-«hllz,2 
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ll cr '- 'hlll,2 


h 3 
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Table 5 
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Test 2b: Change of \\u e - u' h \\ w.r.t. h (£ = 0.001,) 
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